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Abstract 

We have studied the nonhnear current-voltage characteristic of a two dimen- 
sional lattice Coulomb gas by Monte Carlo simulation. We present three 
different determinations of the power-law exponent a(T) of the nonlinear 
current- voltage characteristic, V ~ /'^(^)+^. The determinations rely on both 
equilibrium and non-equilibrium simulations. We find good agreement be- 
tween the different determinations, and our results also agree closely with 
experimental results for Hg-Xe thin film superconductors and for certain sin- 
gle crystal thin-film high temperature superconductors. 
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I. INTRODUCTION 



In two dimensions the superconducting transition in zero magnetic field is a Kosterlitz- 
Thouless transition. This has been verified over the years in both experiments and in 
many models of superconductors like the XY, Villain, and Coulomb gas models • The 

important degrees of freedom in a system undergoing a Kosterlitz-Thouless transition are 
thermally excited vortex pairs. The Kosterlitz-Thouless transition is sometimes also referred 
to as a vortex unbinding transition, as for temperatures below the transition temperature 
Tc all vortices are bound in neutral pairs. These pairs start to unbind at and above T^. 

A typical way to look for a Kosterlitz-Thouless transition in experiments on thin su- 
perconducting films is to probe the current-voltage (IV) characteristic @,|3J^. Both the 
linear and the nonlinear IV characteristics have specific fingerprints identifying a Kosterlitz- 
Thouless transition. Vortices determine the IV characteristic for the following reasons: If 
a vortex is dragged across the system a voltage is induced. Hence resistance is zero only if 
there are no vortices available to move across the system, and only then the system is truly 
superconducting. Vortices that are bound in neutral pairs are unable to move freely and to 
cause dissipation. However an external applied in-plane supercurrent yields a perpendicular 
Lorentz force acting in opposite direction on vortices with different vorticity. This gives a 
net flux of vortices across the system, which shows up as nonlinear (i.e. current dependent) 
resistance. 

Below the Kosterlitz-Thouless transition temperature all vortices are bound in neutral 
pairs by the logarithmic vortex interaction, and the linear resistance is thus zero. Therefore 
the system superconducts below the Kosterlitz-Thouless transition. The linear resistance 
drops to zero at the Kosterlitz-Thouless transition with an exponential functional form, 
R ~ with In^ ~ 1^ ~ Tc]^^^"^ [^. This is consistent with experiments, although the 
logarithm is a complication for quantitative comparison between theory and experiment. A 
finite applied current gives a power-law nonlinear IV characteristic of the form V ~ 
The critical current is thus zero. At the Kosterlitz-Thouless transition the IV exponent a{T) 
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assumes the universal value 2, so ~ at T = Tc. For T < Tc one has a{T) > 2, and for 
T > Tc one has a{T) = (for small enough currents) Experiments on, for example, thin 
Hg-Xe alloy films and also for certain single crystal high temperature superconductors 
I^J^, among some, have confirmed this. 

Since IV characteristics are hard to calculate analytically computer simulation is a useful 
tool. IV characteristics of vortex systems have recently been calculated successfully with 
Monte Carlo simulations |jlO| . Linear and nonlinear IV characteristics of vortex glass super- 
conductors have been reported in Refs. 0,0]. In a recent Monte Carlo simulation of the 



Coulomb gas the linear resistance was used to locate the Kosterlitz-Thouless transition |]I2 
The nonlinear IV characteristics at the Kosterlitz-Thouless transition has been calculated 
in Ref. [jl3|, and a finite-size scaling analysis accurately verified the relation ~ at the 



Kosterlitz-Thouless transition. 

In this paper we study the IV characteristics of a lattice Coulomb gas model by Monte 
Carlo simulations of vortex dynamics. We calculate the IV exponent a(T) of the Coulomb 
gas in three different ways: (1) By direct Monte Carlo calculation of the nonlinear resistance, 
(2) by a self consistent linear screening construction for the energy barrier for current in- 
duced vortex-pair breaking giving thermally activated resistance, and (3) by a finite scaling 
construction from data for the linear resistance. All methods are based on Monte Carlo 
simulations, and we apply both equilibrium and non-equilibrium simulations. These three 
methods give the same results, giving us a consistent and simple picture of nonequilibrium 
response in this system. Furthermore, we compare our results for a{T) with experiments. 
Scaling arguments give that a(T) is a universal scaling function of a reduced Coulomb gas 
temperature X = T/Tc, and this is verified in experiments ||^. We find close agreement 
between our Monte Carlo results and the experimental universal scaling curve, and this 
comparison appears to be presented here for the first time. The agreement between dif- 
ferent methods, and between our simulations and experiments, are the main results of our 
paper. Some of our Monte Carlo results for the nonlinear IV characteristics have been 
obtained previously as explained above. 



The paper is organized as follows: In Section II we define the lattice Coulomb gas 
model. In Section III we study various approaches to the IV characteristics. In Section IV 
we describe our Monte Carlo methods for calculating IV characteristics. In Section V we 
present the Monte Carlo results. Section VI contains discussion and conclusions. 



II. LATTICE COULOMB GAS 

A useful starting point for calculations with superconductors in the presense of currents 
and fields is the Ginsburg-Landau model, with the order parameter \l/(r) = |\I^(r)|e*''^*^'"^ 
describing the superconducting order of the system. However, this model does not focus 
particularly on vortex degrees of freedom. The vortices constitute the essential degrees of 
freedom near the Kosterlitz-Thouless transition. An approximation to the Ginsburg-Landau 
model which focuses only on the vortices is given by the Coulomb gas model. Here thermal 
fluctuations in the magnitude of \I' are neglected, since they are relevant only close to the 
mean-field transition temperature, which is assumed to be well above the vortex transition 
temperature T^. In our simulations the model is discretized and put on a lattice. The 
approximation made in the lattice discretization will only affect the short range behavior of 
the vortices, as the lattice defines the smallest possible separation. The critical properties 
will however not be effected. In general, large length scale properties should be reasonable 
modeled by the lattice Coulomb gas close to Tc. 

The lattice Coulomb gas |1T^JT5[| is defined by the partition function Z on a square lattice 
of side length L using periodic boundary conditions: 

Z = Tr„exp[-/?(ff-/iiV)] (1) 
H = \Y. niGijUj (2) 

i,j 

N = J2\n^l (3) 

i 

where H is the Hamiltonian, Ui is the vorticity at site i (Coulomb gas charge), fi = —Ec is the 
vortex "chemical potential" and Ec is the vortex core energy, and T = is the Coulomb 
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gas temperature 0. The trace is over rii = 0, ±1 on all sites i, subject to overall neutrality, 
J2i = 0. Gij is the lattice Green's function for the logarithmic 2D vortex interaction, 

Gij = J2T. 2 - cos(A;,) - cos{ky) ' ^'^^ 

where k are the reciprocal lattice vectors, k^, ky = 27rn/L, n = 0, . . . , L — 1. 

We will calculate the response voltage to an applied current imposed on the Coulomb 
gas. The above definition does not include any net currents. How to include them and to 
calculate IV characteristics by Monte Carlo simulation is described in the next section. 



III. CURRENT- VOLTAGE CHARACTERISTICS 

In this section we discuss various aspects and approaches to the cur rent- volt age charac- 
teristics of 2D superconductors close to the Kosterhtz-Thouless transition. 



A. Linear resistance 

A basic experiment on a superconductor is to measure the linear resistance. Such mea- 
surements on thin films of both conventional low- Tc superconductors [^] and single crystal 
high-Tc materials have been successfully interpreted in terms of thermally excited 

vortex fluctuations analyzed by use of the Coulomb gas 0. 

The linear resistivity is defined by p = E/j for j — ^ 0, where j is the applied supercurrent 
density and E is the resulting induced electric field. Some words about notation: Since 
resistance and resistivity have the same dimension in two dimensions and our system is 
homogeneous, they are the same, and they we will both be denoted by R. R will be 
reserved for linear resistance, and will not be used to denote nonlinear resistance. An 
applied supercurrent is denoted by J = jL, and voltage is ^ = EL. 

To determine the linear resistance in simulations of the Coulomb gas from E/j for small 
j has its limitations, as we have to repeat the calculation at a sequence of current densities 
j, to make sure that j is small enough to be in the linear regime. If the purpose is to 
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measure only the linear resistance, and not E as function of j, a different approach is to 



use the Nyquist formula [0, which relates the linear resistance to the equilibrium voltage 
fluctuations: 

1 f+oo 

where V{t) is the induced voltage from vortex motion at time t. As an alternative to eq. 
H the Kubo formula for the vortex currents 1^, R = ^ I^oo^'^ {Iv{t)Iv{0)) can be used. 
Given the Josephson relation we see immediately that the Kubo formula equals the Nyquist 
relation. 



The linear resistance has been successfully used in a simulation [|T2[ to locate the 

Kosterlitz-Thouless transition temperature Tc of the 2D lattice Coulomb gas. They find 
the finite size scaling relation at T^. 

L^R ( 1 + J = constant at T = T^, (6) 

to be valid to a very high precision. 

The scaling relation Eq. was derived from the following argument. We assume the 

dynamical exponent z = 2 for free vortex diffusion in two dimensions []T^|T^. The linear 



resistance is a dynamical quantity, it relates to the correlation time r, which at Tc diverges 
like r ~ where ^ is the correlation length. According to the Josephson relation the 
voltage V ~ ~ r~^, where A(f) is the gradient of the phase of the Ginsburg-Landau 
order parameter |jl8[. Therefore, we expect the linear resistance, Eq. to scale like 



R ~ ^"^ at Tc. At Tc the correlation length diverges and is cut of by the finite size L of the 
lattice and hence RL^ = const at Tc, to lowest order. The scaling relation has a logarithmic 
correction which has been included in Eq. (^. This correction is readily obtained from the 
corresponding correction terms for 1/e and A ||T^ . 
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B. Thermally activated resistance 



The above scaling argument led to a finite size scaling formula which is useful for locating 
the transition temperature from Monte Carlo data for the resistance of finite samples. Here 
we will do a more detailed analysis that will also lead to the same formula. The analysis here 
does not directly involve scaling arguments, but considers the interactions between vortices 
in the Coulomb gas. The analysis will give expressions for the resistance from thermally 
activated free vortices in the Coulomb gas in the presence of an applied supercurrent. This 
more detailed analysis will be useful in later sections when we analyze Monte Carlo data for 
the Coulomb gas. 

According to the Josephson relation the voltage V caused by vortex motion is 

V ~ ~ npl 

dt 

where we assume that the resistance is proportional to the density of "free" vortices, np, 
defined by the the Debye-Hiickel relation. The linear resistance R is defined by the limit of 
zero current J: 

V 

R = lim — Up (7) 

i~*o I ^ ' 

To make an estimate of the density of free vortices we proceed by the following simple model. 



The energy E{r) of a vortex pair of separation r > rg in the presence of a current / is |T7 



E{r) =Eo + E^ In(-) - /(r - tq) (8) 
f 

where Eq is a constant, Ei is discussed below and Tq is the smallest possible separation, 
which we will set to tq = 1 from now. 

We will now use the linear screening approximation to derive an expression for the 
second term Eilnr in Eq. (^. The expression is obtained from the Fourier transform of 
the linearly screened potential Vi{k), 

1 2ti 



Vi{k) 



i{k) P + A-2 
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Here A is the vortex screening length, and e{k) is the part of the dielectric function, e{k), 
describing the polarization of the bound pairs. The two e are related by 

1 1 P 



e(k) e(k) P + A-2 

In the limit A — oo the two different dielectric functions become equal. This is the case 
for temperatures below T^. The dielectric function, l/e{k), is obtained from the charge 
fluctuations below in Eq. (p5l). The real space expression for E{r) is obtained from 

E{r) = lim Vi{r) 

A^oo 

Where 



We can obtain an approximate expression for Vi{r) by making use of the fact that e(k) only 
depends weakly (in most of k space) on k. For a given distance r, the Fourier integral picks 
up its main contribution from the k values around 27r/r. Hence 

V1(r) - mr = 1) « -^^^J^^K,(r/Xy 

Here we have subtracted Vi{r = 1) in order to eliminate the creation energy. Kq denotes a 
modified Bessel function. As A ^ oo this expression reduces to 

V^(r)-V^(r = l)^— 1— ln(r/A) (9) 
e(27r/r) 

where we use e instead of e, as the temperature is below Tc. 
According to this discussion the coefficient Ei is given by pO 



El = , \ (10) 

e(27r/r) ^ ^ 

The weak r dependence describes the effect of the surrounding vortex pairs. The coefficient 
Eq contains the remaining constant terms from Eq. (^. In a first approximation we will 
neglect the r dependence in Ei. 
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The energy E{r) in Eq. (§) has a maximum at separation r* = Ei/I and the energy 
needed to separate a vortex pair to this distance is |^ : 

AE = E{r*) - E{r = I) = Ei ln(r*) - /(r* - 1) (11) 

AE = Ei\n{^)-Ei + I (12) 

Let r denote the thermal production rate of free vortices. A vortex vanishes when it 
colhdes with an antivortex. Hence, it appears reasonable to assume an annihilation rate 
proportional to n^. This leads to the following rate equation. 

fij? = r — cn?p. (13) 



Where c is a constant. The steady state condition is hp = ^ and hence we have F oc y/np. 
Assuming that F is determined by activation over the barrier AE we get the following 



production rate [|T7 



F oc e~— (14) 
and hence for the resistance R from Eq. (|7]). 



AS 



R (X Up OC e 2T oc 

.El.Ei £i i_ 

^ j 2T g2Tg 2T 

keeping the important term for small but finite / we arrive at: 

R(x{j-)'^ (15) 

A given current / gives rise to a "current length scale" r* from the maximum condition in 
Eq. (H). As the lattice of the system has a finite size, this sets an upper limit to the "current 
length" and hence a lower limit to the current producing nonlinear resistance. The smallest 
current giving nonlinear resistance is /* = Ei/r* with r* = L and hence for currents smaller 
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than /* the resistance will be cut off by the finite size L of the lattice and the resistance 
becomes ohmic. The Nyquist resistance is calculated with 7 = and hence 

R(x{]r)-^. (16) 

This means that we can scale the linear resistance R from the Nyquist relation Eq. 
with the exponent Ei/2T. This exponent is precisely a(T), the exponent of the nonlinear 
IV characteristics (see Eq. (|TEp below), hence: 

/(T) = i?L"(^) (17) 

should collapse onto a single curve for different lattice sizes L. I.e. /(T) should not depend 
on lattice size L. The resistance we use for this scaling will be the one determined from the 
voltage fluctuations Eq. (^. The exponent determined from resistance data at zero current 
will be denoted a/j(T). 

C. Nonlinear IV exponent 

We are going to make use of a couple of different expressions for the power law expo- 
nent a(T) of the nonlinear IV characteristics. From Eq. (|T5|) we get the nonlinear IV 
characteristic: 

y oc (^)-i^Joc (18) 

The exponent calculated by monitoring the voltage response \^ as a function of an applied 
supercurrent / will be denoted a/y(T). On a finite system we will obtain a nonlinear voltage 
response only above a finite applied current, given by /* ~ Ei/L, such that the current length 
r* is shorter than the size L of the system, as discussed above. 

D. Self consistent IV characteristic 

Another expression for the IV characteristic is obtained if we include the r dependence 
in El in Eq. ([lO|) . The length dependence can in a first approximation (in an expansion in 
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derivatives of Ei{r)) be included simply by replacing Ei in Eq. (|T5|) by l/e(27r/r*) in the 
extremum equation / = Ei/r*. Our rationale for this choice is that at the separation r* the 
vortex pair is broken apart and we therefore use the stiffness l/e(r) of the system at this 
separation. We find the appropriate e(r) by solving self consistently the equation 

1 k* 

^ ~ e(27r/r*)r* ~ e(A;*)27r ^^^"^ 

The self consistent e obtained by solving Eq. ([ToD will be denoted e*. The relation between 
the exponent a(T) and the dielectric function e is according to Eqs. (|15|) and ([T0|) given by 
the expression (see Ambegoakar et al. [p!7| ) 

a{T)AHNS = ^ (20) 

here we use e* as we are at temperatures below Tc. 

Recently Minnhagen et al. have used scaling arguments to derive an alternative relation- 
ship between a(T) and e, given by ||21[ 



a{T)pM = ^ - 2 (21) 

As one immediately realises Eq. (^) is not consistent with the activation argument used 
to derive Eq. (pO]). In order to reconcile Eq. ( ^I]) with a rate equation like Eq. ([13D 



Minnhagen et al. have made the following suggestion. They assume that the activation is 
correctly represented by F in Eq. ([T^). The recombination, which in Eq. ( p!3D is represented 
by the innocently looking term n|,, is on the other hand supposed to be replaced by n]i^'^ 
with b = 2/{Ei/T — 2). The sole argument for this replacement is unfortunately so far 
simply the observation that one then can derive Eq. (pT]) from an equation like Eq. (|13|). 



Nonetheless, we shall see below that for temperatures below Eq. (|21| ) fits the simulation 
data much better than Eq. ( pOf ) does. However a motivation for a recombination term 
different from the one in Eq. (|13[) has not been presented. At Tc both relations reproduce 
the same exponent a{T = Tc) = 2. 
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IV. MONTE CARLO SIMULATION 



In this section we describe how we calculate current-voltage characteristics by Monte 
Carlo simulation of the lattice Coulomb gas. 



The algorithm to simulate the lattice Coulomb gas works as follows [15|: First we pick 
a nearest-neighbor pair of lattice sites at random. Then we try to increase Ui by one 
and to decrease rij by one, thus preserving overall vortex neutrality, X^i^i = 0. This Monte 
Carlo move of inserting a neutral pair will be interpreted as transfer of one unit vortex from 
site j to i. If the energy change is AE we accept this trial move according to the standard 



Metropolis algorithm [^] with probability exp{—AE/T). These simple Monte Carlo moves 
can both create, annihilate, and move vortices. Thermodynamic averages are computed as 
Monte Carlo time averages over the sequence of generated configurations. 



To calculate IV characteristics works as follows |TO|JTl|: An applied current density j 
gives a Lorentz force of jh/{2e) on a unit vortex. The Lorentz force can be incorporated in 
the Monte Carlo moves [jlO[ by adding to AE an extra term jh/{2e) if the unit vortex moves 



in the direction opposite to the Lorentz force, subtracting this term if it moves in the same 
direction, and making no change in AE if it moves in a perpendicular direction. Biasing the 
Monte Carlo moves in this way takes the system out of equilibrium and causes a net flux of 
vortices in a direction perpendicular to the current. This generates a voltage given by the 
Josephson relation: 

V = ^jm), (22) 

where /^(t) is the vortex current. Here t denotes Monte Carlo time, incremented by 6t after 
each attempted move. The vortex current is Iv{t) = +l/LAt if a unit vortex moves one 
lattice spacing in the direction of the Lorentz force at time t, and Iv{t) = otherwise. We 
use units such that At = l/L"^ so that an attempt is made to move a vortex on each lattice 
site, on average, per unit time. We also use units such that h/{2e) = 1. 

The linear resistance can also be obtained from the Nyquist relation in Eq. (13) for 
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equilibrium voltage fluctuations in the absence of any net currents. For discrete Monte 



Carlo time it is given by |T6,23|, in our units, 



^ = ^ E Mvmrn. (23) 

t=-T 

The cutoff time r is set to 1000 time steps, this has proved to be sufficient as R saturates 
for values of r less than 100 time steps. 

The underlying assumption in using Monte Carlo dynamics to calculate IV characteristics 
is that Monte Carlo time can be equated with real time. This approximation has proven 
reasonable in other simulations of vortex dynamics pO| jTl|JT3[1 . It should be good near a 
critical point where vortex motion is slow and overdamped, but not so satisfactory at high 
temperatures or currents where the discreteness of Monte Carlo time becomes visible as a 
saturation of vortex velocities. We get strong support for this assumption from the results 
in the next section since we can reproduce the expected IV characteristic at the Kosterlitz- 
Thouless transition, and since we come close to experiments. One can in principle also test 
this by comparison to dynamics simulations where time evolution equations are integrated 

g. 

In the equilibrium simulations in the case of no net currents we typically use 10^ — 10^ 
Monte Carlo sweeps (one sweep means one Monte Carlo time step deffned above, i.e. L x L 
attempts to insert nearest-neighbor pairs), and in the nonequilibrium case of an applied 
current we typically use 10^ — 10^ Monte Carlo sweeps. In the evaluation of the Nyquist 
formula for the linear resistance, Eq. (^, we typically sum over 10^ — 10^ time steps. 

The ffrst task for our simulations is to locate the Kosterlitz-Thouless transition temper- 
ature Tc. The usual universal jump criterion for a Kosterlitz-Thouless transition involves 

the dielectric response function 1/e, given by 

1 2tt 

= (25) 

where nk is the Fourier transform of the vortex density. The limit k — denoted l/e(k = 0), 
corresponds to the fully renormalized long wavelength superffuid density, and the universal 
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jump criterion tells us that l/e(k = 0) jumps from 4Tc at T = T~ to at T = T+ 0. 
A practical difficulty for locating from Monte Carlo data on small lattices with this 
procedure is that extrapolation to the k = limit requires large lattices, as the smallest 
nonzero k is 2tt/L. The corresponding quantity to l/e(k = 0) in the two dimensional XY 
model is called the helicity modulus 7 0. Both quantities have been used to locate the 
Kosterlitz-Thouless transition temperature in Monte Carlo calculations [|5|, p5| , [T9 |. 

In the data analysis in the next section we use an alternative procedure to locate Tc from 
the linear resistance [|12|. We obtain the linear resistance R from the Nyquist formula in 
Eq. for a sequence of system sizes L and temperatures T. According to Eq. @ data 
for L^R for different system sizes should become system size independent at the critical 
temperature, which is our criterion to locate T^. This circumvents the difficulties of using 
l/e(k = 0). The two determinations give within error bars the same value for T^. 



V. RESULTS 

We present Monte Carlo simulation results for three different determinations of the IV 
exponent a{T) for the 2D Coulomb gas model. The results we show here are for the chemical 
potential fi = 0.0 and lattice size L = 32 if not differently stated. 



A. Nonlinear IV exponent 

Our ffist method consists in a direct measurement of the electric field E induced by an 
applied current density j. In Fig. |l]a results for the IV characteristic of the two-dimensional 
Coulomb gas are shown. The dashed line in the ln{E) versus ln(j) plot has slope three 
and represents the slope at T = according to the universal jump condition |p. The 
solid curves represent results for different temperatures. For very high current the voltage 
response saturates. This is because when all attempts to move the vortices in the direction 
of the Lorentz force are already accepted, further increasing the current can not give more 
voltage. For low enough current there is a crossover to ohmic resistance, when the current 
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length equals the system size, and the nonlinear dependence of the resistance on the current 
vanishes. The regime where we probe the non-linear IV characteristic is for this figure 
approximately from Inj ^ —1.5 up to ^ —0.5. According to the Kosterlitz-Thouless theory 
the slope of the lines should be 3 at the critical temperature, and this criteria can be used to 
determine Tc. We will however use an independent determination of for this system, 
based on the finite size scaling relation Eq. (^). 

In Fig. |l]b we demonstrate the effects of the finite lattice size for low driving currents 
at temperatures below Tc. The data shown are for T = 0.15 and lattice sizes are L = 8 
(triangles), 12 (open squares), 16 (stars), 24 (open circles), and 32 (filled circles). 

The finite size effects for the lower temperatures can be understood in the following way. 



(See the discussion above following Eq. (pIS]).) The finite lattice size is important because 
pair excitation over the barrier given by the periodicity length L will add to the dissipation 
due to unbinding of pairs over the barrier given by the pair size r*. The induced electric 
field will accordingly be of the form 

E = R{L)j + constant (26) 

where the first term R{L) follows from Eq. ( [T^ ) and R{L) — as L — oo. The second 
term in Eq. (|26|) is given by Eq. ([ISD and will remain finite in the limit L oo. This is 
clearly demonstrated in Fig. 03 where we see that the crossover in Eq. (|26|) between the 
linear and nonlinear regime appears at a higher driving current for the smaller 8x8 lattice 
as the current length {Ei/I = ^/ ~ L) associated with the current density j exceeds the size 
of the lattice. 

In Fig. the exponent a/y(X) is shown as a function of the reduced temperature 
X = T/Tc. The dashed horizontal line represents the universal jump condition for ajv{X). 
The plusses represent experimental data from a superconducting Hg-Xe film P,p6|. The 
filled circles are the results for aiv{T) from Fig. |1]. The other three data sets are for 
lattice sizes L = 16 (stars), 24 (open circles), and 48 (triangles). As one can see there are 
no apparent finite size effects in the data. In the vicinity of the critical temperature the 
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experimental data are reproduced by the Monte Carlo simulations. 

The reduced temperature variable X used in Fig. |^ for the experiment is also from Ref. 
p6| and for the lattice Coulomb gas data we use Tc = 0.218 [|T^ determined from a finite 
size scaling analysis using Eq. (||). 

The inset in Fig. ^ shows a selection of experimental data analyzed along the lines 



described in Ref. [^. The data in the inset (plusses) are the same as in the main 
figure, the other data are for Bi2Sr2CaCu20a; single crystal (filled squares) 0, and for 
Bii.6Pbo.4Sr2Ca2Cu302^ single crystal (open squares) p|. 

The Monte Carlo data presented here for a{X) are all for /i = 0.0. We also did the same 
analysis for Monte Carlo data for L = 32 and /i = —0.4, —0.2 and 0.2. The closest fit to the 
experimental results is produced hj ji = 0.0. Results for different /i differ from the n = 0.0 
results, by that /i = 0.2 has a slightly larger derivative at X = 1 and the smaller /i are 
correspondingly less steep. 

B. Self consistent IV characteristic 

In our second determination of the exponent a{T) we will make use of the relations 
between e* and a{T) in equations (|20| - |2T1) . The analysis is based on the self consistent 



solution of Eq. (0). For a given current density j, a set of e{k) will be calculated for 
different temperatures. The self consistent solution to Eq. (|T9|) for e* is shown in Fig. |^, 
in (a) data for T = 0.18 is shown and in (b) T = 0.24. The solid and dashed straight lines 
represent e* = fc*/j27r, given by Eq. (|l^), for different current densities j. The open circles 
represent e(k*) as a function of k for different current densities. The choice of the direction 
along which e(k*) is probed is perpendicular to the current density j, ie. parallel to the 
vortex drift caused by the current density j. The intersection between a e(k*) curve and the 
corresponding straight line is the solution to Eq. (0), these are marked with filled circles. 
In Fig. ^ for T = 0.18 < Tc = 0.218 we see that the self consistent solution depends only 
weakly on the choice of probing current as long as the current is not too large. In Fig. ^ 
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for T = 0.24 > Tc we see, however, that there is no well defined limiting solution for e* as 
j ^ 0. This is because the system is above the Kosterlitz-Thouless temperature and vortex 
pairs will always dissociate irrespective of j. In Fig. ^ the function e* is shown as a function 
of temperature. The solid circles represent e* from the self consistent Eq. (|T9D for the fixed 
current density j = 0.03125. The data shown here represents the construction shown in Fig. 

i 

The results from the self consistent solution for e* are analyzed in Fig. |^. Here the filled 
circles represent the exponent a/y(T) from Fig. |^. The upside down triangles represent 
aAHNs{T) from Eq. (|20D with the solution from Fig. ^ and the triangles are the corre- 
sponding solution to Eq. (^). One can clearly see that the expression in Eq. (0), derived 
by Minnhagen et al. [^, reproduces the exponent ajv{T) for T < T^. Note however, it is 
only a coincidence that Eq. (pO]), derived by Ambegoakar et al. |1^, works for temperatures 
above Tc in this figure as the limiting (j — > 0) solution for 1/e* is not well defined for these 
temperatures, as already discussed in connection with Fig. |]b. As the simulation data 
aiv{T) (filled circles) also matched the experimental data in Fig. |I] we must conclude that 
below Tc the interpretation according to Eq. (^TJ) is clearly the more appropriate. 



C. Linear Resistance 

We will now turn to our last determination of the exponent a(T). The results presented 
above all relied on non-equilibrium Monte Carlo simulations, i.e. with a finite applied 
supercurrent density j. We will now present the equilibrium determination for j = based on 
finite size scaling of Monte Carlo data for the linear resistance given by the Nyquist formula 
(I) together with Eq. ([TtD . In Fig. ^jwe demonstrate a data collapse of the hnear resistance 
for several lattice sizes. From Eq. (p!7| ) we see that the linear resistance data can be 
collapsed onto a single curve, thus representing the thermodynamic limit, by an appropriate 
choice at each temperature T of the exponent a/j(T). We do this in the following way. For 
a given temperature we find the exponent which minimizes the error of the fit defined 
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as EL^L'iRWL" - R{L')L'^f. The considered lattice sizes are L,L' = 6,8,12,16,24,32. 
The obtained scahng exponent as a function of temperature is shown in Fig. ^ As 
a comparison we also show data for a/y(T) from Fig. |^, obtained from direct evaluation 
of the IV characteristic. The finite size scaling analysis in Fig. ^ breaks down for low 
temperatures. This can be seen by the deviation of aji{T) from the data for ajv{T) at 
T = 0.15. In Fig. this deviat ion is also evident. A careful inspection of the scaling at 
temperatures T = 0.15 and T = 0.18 reveals that the order of the lattices sizes is reversed 
for T = 0.15 compared with the higher temperatures. This may be related to the difficulties 
to converge the simulation at low temperature. 



VI. DISCUSSION 

We have calculated the nonlinear IV exponent ajv{T) of the two dimensional lattice 
Coulomb gas. Our results are based on three different determinations. A direct calculation 
of the voltage response as a function of an applied current. Comparison with experiments 
P|7|J^,P^ on Hg-Xe films and single crystal high superconductors show good agreement. 

Our second method is based on a simple self consistent calculation of the dielectric 
function e* at the unbinding separation, and the IV exponent can then be calculated. Here 
we especially focus on the comparison of two relations between a(T) and e. The first relation 
Eq. (^) |]l^ is based on ordinary diffusion in two dimensions with a recombination rate 
proportional to rPp. The second expression for a(T) given in Eq. (^) has been derived from 
a scaling analysis ||21[ . 



We find that the exponent determined by Eq. (^) for temperatures below is close 
to the more direct determined aiv{T) and will therefore also fit the experiments for these 
temperatures. 

The third method is based on equilibrium Monte Carlo simulations. From the scaling 
relation Eq. ( [T7| ) for the linear resistance we can derive afi{T). We find that the scaling 
exponent aij(T) to a high degree of accuracy fits the direct determined a/y(T) for a broad 
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range of temperatures. This provides a link between the equihbrium and nonequihbrium 
response properties of the system: A finite mesoscopic hnear resistance in a finite sample 
below Tc is due to thermally activated vortex motion across some potential barrier, generated 
by interactions with all other vortices in the system. When a finite current is imposed across 
the system, new nonequihbrium configurations are accessed where vortices are driven away 
from their equilibrium positions by the finite Lorentz force, thus giving nonlinear response. 
Our data shows that the barrier overcome by the Lorentz force, giving nonlinear response, 
is essentially the same barrier as in the equilibrium case, i.e., the potential barrier in the 
case of a finite current appears to be determined by equilibrium states in the system. This 
is consistent with the scaling ansatz, discussed above, of a certain current length scale (r*) 
associated with the finite current, such that for lengths shorter than the current length scale 
an equilibrium state is still attained, which gives a potential barrier essentially equal to that 
in the equilibrium case. 

An interesting possibility arises here to measure finite size effects on the linear resistance 
in lithographic Josephson junction arrays. The idea would here to take advantage of finite 
size roundings, rather than as usual want them to be as small as possible, and trying to 
fit experimental data on very small arrays to our finite size scaling formulas. This would 
provide an unusual experimental test of finite size scaling. It is also important to analyze 
the data in terms of the reduced temperature scale, as Tc is sample dependent. According to 
Eq. {^n\) it should be possible to scale the linear resistance of samples of different sizes onto 
a single scaling function using the exponent a{X) from the nonlinear IV characteristics. 

Our conclusions are: (1) Simulation of nonequihbrium vortex dynamics allow calculation 
of a lattice size independent IV exponent a(T) as a function of temperature T. (2) This 
curve for a(T) agrees nicely with experiments in an interval around T^, and this appears to 
be reported here for the first time. (3) This curve can be obtained from a simple phenomeno- 
logical theory for the nonlinear IV characteristic. (4) This curve can also be obtained from 
a simulation of the equilibrium vortex dynamics. This provides a useful link between driven 
diffusion and equilibrium dynamics of two-dimensional vortex systems. 
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FIGURES 

FIG. 1. Monte Carlo results for the nonlinear IV characteristic of the two dimensional lattice 
Coulomb gas. E is the electrical field and j is the supercurrent density. In (a) data shown are 
for parameters L = 32 and /i = 0.0. Different curves are for different temperatures, starting 
from below T = 0.12, .15, .16, .18, .20, .22, .23, .24, .26, and .30. The dashed line has slope 3 which 
represents the Kosterlitz-Thouless transition. In (b) finite size effects according to Eq. ( |26|) are 
demonstrated. Data here is for T = 0.15 and lattice sizes are L = 8 (triangles), 12 (open squares), 
16 (stars), 24 (open circles), and 32 (filled circles). The dashed line has slope 1 which represents 
the linear ohmic regime. 

FIG. 2. Comparison of the exponent a{X) between experiments, marked with plusses p6| , 
and Monte Carlo results for the lattice Coulomb gas Sizes are L = 16 (stars), 24 (open circles), 
32 (filled circles), and 48 (triangles). X = T/Tc is the reduced Coulomb gas temperature. The 
dashed line corresponds to the universal jump condition for the exponent, a{X) = 2. The inset 
shows experimental data for Hg-Xe alloy films (plusses) [§J2^, Bi2Sr2CaCu20a; single crystal films 
(filled squares) [^, and Bii,6Pbo.4Sr2Ca2Cu30x single crystal films (open squares) [^. 

FIG. 3. Construction of the self consistent solution. The data shown is for T = 0.18 in (a) 
and T = 0.240 in (b). The different curves represent e(/c) as a function of k^. for different current 
densities j = 0.05 (bottom curve), 0.10,0.15,0.20,0.30 (top curve). The straight lines represent the 
self-consistency condition for the different current densities, and the intersections of the straight 
lines with the €{k) curve for the same current density j represents the self consistent solution, here 
marked with large filled circles. 

FIG. 4. The solution to the self consistent equation (|l9|). e* as a function of temperature for 
the fixed current density j = 0.03125. Data shown are for size L = 32 chemical potential = 0.0. 
The dashed line represents the bare e = 1. 
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FIG. 5. The IV exponent a{T) as a function of temperature. Data shown are for size L = 32 
and chemical potential = 0.0. The filled circles are ajv(T) from the the nonlinear IV character- 
istic shown in Fig. [l|b. The upside down triangles is the exponent aAHNs{T) from Eq. (pO|) using 
e* from Fig. (Q). The triangles represent the exponent apM{T) from Eq. (^T|) also using e* from 
Fig. |. 

FIG. 6. Finite size scaling of the linear resistance R. The resistance has been calculated with 
the Nyquist formula (|5|). Shown is a data collapse of i?L'"«(^) as a function of T for different lattice 
sizes, L = 6,8,12,16,24,32. The exponent aji{T) is chosen in such a way as to provide the best 
data collapse. According to Eq. ([T7|) the function RL""^^^^ should be independent of system size. 

FIG. 7. Comparison of the exponent a{T) as a function of temperature T, from two different 
determinations. The filled circles represent the exponent aiv{T) determined in Fig. ^ from the 
nonlinear IV characteristics. The open circles are results for aij(T) from Fig. |6| from finite size 
scaling of the linear resistance. 
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